/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of LINAL (LINear ALgebra)
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include "linal.hxx"
#include "flops.h"
#define LINAL_ERROR_TOL 1.0e-12

#define LINAL_MUL(is_complex) ( (is_complex) ?  (6.0) : (1.0) )
#define LINAL_ADD(is_complex) ( (is_complex) ?  (2.0) : (1.0) )
#define LINAL_MEM(is_complex) ( (is_complex) ? (16.0) : (8.0) )


int main(int argc, char **argv) {

  FLA_Init();
  
  double flop_linal, flop_magma;
  
  for (int is_complex=0;is_complex<2;++is_complex) {
    for (int k3=0;k3<3;++k3) {
      for (int k2=0;k2<3;++k2) {
        for (int k1=0;k1<3;++k1) {
          double dk1 = (double)k1;
          double dk2 = (double)k2;
          double dk3 = (double)k3;

          // gemm
          flop_linal = linal::get_flop_gemm(is_complex, k1, k2, k3);
          flop_magma = ( LINAL_MUL(is_complex)*FMULS_GEMM(dk1,dk2,dk3) +
                         LINAL_ADD(is_complex)*FADDS_GEMM(dk1,dk2,dk3) );

          if (abs(flop_linal-flop_magma) > LINAL_ERROR_TOL )
            printf("gemm %d %d %d\n", k1,k2,k3);
          
          // syrk
          flop_linal = linal::get_flop_syrk(is_complex, k1, k2);
          flop_magma = ( LINAL_MUL(is_complex)*FMULS_SYRK(dk1,dk2) +
                         LINAL_ADD(is_complex)*FADDS_SYRK(dk1,dk2) );

          if ( abs(flop_linal-flop_magma) > LINAL_ERROR_TOL )
            printf("syrk %d %d %d\n", k1,k2,k3);

          // trsm_lower
          flop_linal = linal::get_flop_trsm_lower(is_complex, k1, k2);
          flop_magma = ( LINAL_MUL(is_complex)*FMULS_TRSM_L(dk1,dk2) +
                         LINAL_ADD(is_complex)*FADDS_TRSM_L(dk1,dk2) );

          if ( abs(flop_linal-flop_magma) > LINAL_ERROR_TOL )
            printf("trsm lower %d %d %d\n", k1,k2,k3);

          // trsm_upper
          flop_linal = linal::get_flop_trsm_upper(is_complex, k1, k2);
          flop_magma = ( LINAL_MUL(is_complex)*FMULS_TRSM_R(dk1,dk2) +
                         LINAL_ADD(is_complex)*FADDS_TRSM_R(dk1,dk2) );

          if ( abs(flop_linal-flop_magma) > LINAL_ERROR_TOL )
            printf("trsm upper %d %d %d\n", k1,k2,k3);

          // trmm_lower
          flop_linal = linal::get_flop_trmm_lower(is_complex, k1, k2);
          flop_magma = ( LINAL_MUL(is_complex)*FMULS_TRMM_L(dk1,dk2) +
                         LINAL_ADD(is_complex)*FADDS_TRMM_L(dk1,dk2) );

          if ( abs(flop_linal-flop_magma) > LINAL_ERROR_TOL )
            printf("trmm lower %d %d %d\n", k1,k2,k3);

          // trmm_upper
          flop_linal = linal::get_flop_trmm_upper(is_complex, k1, k2);
          flop_magma = ( LINAL_MUL(is_complex)*FMULS_TRMM_R(dk1,dk2) +
                         LINAL_ADD(is_complex)*FADDS_TRMM_R(dk1,dk2) );

          if ( abs(flop_linal-flop_magma) > LINAL_ERROR_TOL )
            printf("trmm upper %d %d %d\n", k1,k2,k3);

          // lu
          flop_linal = linal::get_flop_lu(is_complex, k1, k2);
          flop_magma = ( LINAL_MUL(is_complex)*FMULS_GETRF(dk1,dk2) +
                         LINAL_ADD(is_complex)*FADDS_GETRF(dk1,dk2) );

          if ( abs(flop_linal-flop_magma) > LINAL_ERROR_TOL ) 
            printf("%d lu %d %d %d\n", is_complex, k1,k2,k3);


          // chol
          flop_linal = linal::get_flop_chol(is_complex, k1);
          flop_magma = ( LINAL_MUL(is_complex)*FMULS_POTRF(dk1) +
                         LINAL_ADD(is_complex)*FADDS_POTRF(dk1) );

          if ( abs(flop_linal-flop_magma) > LINAL_ERROR_TOL )
            printf("chol %d %d %d\n", k1,k2,k3);

          // QR
          flop_linal = linal::get_flop_qr(is_complex, k1, k2);
          flop_magma = ( LINAL_MUL(is_complex)*FMULS_GEQRF(dk1,dk2) +
                         LINAL_ADD(is_complex)*FADDS_GEQRF(dk1,dk2) );

          if ( abs(flop_linal-flop_magma) > LINAL_ERROR_TOL )
            printf("qr %d %d %d\n", k1,k2,k3);

          // Generating Q
          flop_linal = linal::get_flop_q(is_complex, k1, k2, k3);
          flop_magma = ( LINAL_MUL(is_complex)*FMULS_UNGQR(dk1,dk2,dk3) +
                         LINAL_ADD(is_complex)*FADDS_UNGQR(dk1,dk2,dk3) );

          if ( abs(flop_linal-flop_magma) > LINAL_ERROR_TOL )
            printf("q %d %d %d\n", k1,k2,k3);

        }
      }
    }
  }

  FLA_Finalize();

  return -1;
}
